Cascaded / Unary union #1246

urschrei · 2024-10-31T16:43:14Z

I agree to follow the project's code of conduct.
I added an entry to CHANGES.md if knowledge of this change could be valuable to users.

~~This is a draft since it currently has unknown perf characteristics, and there may be better ways of defining the trait implementers~~

lnicola · 2024-10-31T17:24:43Z

Did you ask for some polygons? They're not WKT (I can help with the conversion to if needed), but here's a bunch of them I just happened to have on hand: https://data.europa.eu/data/datasets/d18f6f95-1d96-4b86-a4ba-291a416e45bb?locale=en.

They're a mix of polygons and multipolygons, but I'm not sure why that would be a problem.

geo/src/algorithm/bool_ops/mod.rs

urschrei · 2024-10-31T22:10:59Z

Some preliminary perf data:

a 73.8 mb WKT file, containing 39490 polygons (extract from large-scale agricultural data, thanks @lnicola!) can be unioned into a MultiPolygon in around 11 seconds in release mode. The naïve approach:

let unioned = data
    .iter()
    .fold(MultiPolygon::<f64>::new(vec![]), |acc, poly| {
        acc.union(poly)
    });

hadn't terminated after 18 minutes, so this is at least a 100x improvement (?).

Next I need to verify that the algorithm produces valid output (no reason to think it doesn't)

cc @dabreegster please take this for a spin with more real-world data when you get a chance!

urschrei · 2024-11-01T11:14:18Z

Since we can union MultiPolygons, I'd like to make this available for them too, but I don't think it's possible to generically express that, because you need concrete types for the init, fold, and reduce function definitions and they vary slightly between Polygon and MultiPolygon.

I'd also be happy to write a second impl for MultiPolygon, but rustc also complains about that for some reason, even though the bounds are different:

impl<T, Container> UnaryUnion for Container
where
    T: BoolOpsNum,
    Container: IntoIterator<Item = MultiPolygon<T>> + Clone,
    MultiPolygon<T>: RTreeObject,
{ ... }

dabreegster · 2024-11-01T17:20:44Z

I'm seeing great results from INSPIRE parcels. https://www.dropbox.com/scl/fi/8w9ufilq2c7wq1vp4fei0/London_Borough_of_Southwark.zip?rlkey=ejcug1i5fnkf35d2xjwcogz9m&dl=0

Using mapshaper v1.geojson -dissolve -explode -o v2.geojson format=geojson geojson-type=FeatureCollection, it takes 1.3s to dissolve 40k features into 7.5k. Using the new API in this PR (via https://github.com/acteng/will-it-fit/tree/use_unary_union), it takes 5.1s to dissolve to 1.1k. This implementation handles inputs that already overlap better than mapshaper.

Mapshaper output, darker area is features overlapping originally that weren't dissolved:

Rust output handles these:

There's more, larger INSPIRE inputs. I'm on battery power right now, so not going to push it :) But a larger run in Birmingham on 60k features ran in about 30s (I didn't time it exactly) and produced 17k outputs. I have even larger OS data I could try, but can't share the inputs.
https://github.com/acteng/will-it-fit/blob/use_unary_union/data_prep/inspire_one_area.sh

It's particularly impressive/useful that the base boolean ops work on both projected and un-projected inputs!

No feedback on this API; it seems fine to me. I have a more niche use case of limiting the maximum area of unioned polygons, so I'll probably adapt what's here to do that. I don't think baking in support for that level of customization makes sense.

michaelkirk · 2024-11-01T17:32:20Z

un-projected inputs

As I understand it, the math is all fundamentally euclidean, so you're playing with fire if you're doing these calculations without first projecting to the euclidean plane.

Fire can be kind of fun though.

michaelkirk · 2024-11-01T17:34:29Z

speaking of a multithreading flag - this feature might be a good candidate for rayon. (probably not in this PR, but it'd be interesting to see)

kylebarron · 2024-11-01T18:20:24Z

this feature might be a good candidate for rayon

It might be good to have a separate discussion thread about rayon API design. In particular, I'd love it if any rayon integration could use a rayon thread pool provided by the caller, so that higher level apps that already use rayon don't have work split across multiple rayon pools.

michaelkirk · 2024-11-01T23:27:49Z

Since we can union MultiPolygons, I'd like to make this available for them too,

Could you implement it in terms of impl BooleanOps<Scalar = Self::Scalar> (both MultiPolygon and Polygon implement that)

Something like:

- Container: IntoIterator<Item = Polygon<T>> + Clone,
+ Container: IntoIterator<Item = impl BooleanOps<Scalar = T>> + Clone,

This should work for both Iterator<Item=Polygon> and Iterator<Item=MultiPolygon>, but if you wanted a single iterator, you'd need to do something like Iterator<Item=Box<dyn BooleanOps>>

urschrei · 2024-11-02T00:07:18Z

Since we can union MultiPolygons, I'd like to make this available for them too,

Could you implement it in terms of impl BooleanOps<Scalar = Self::Scalar> (both MultiPolygon and Polygon implement that)

Something like:
- Container: IntoIterator<Item = Polygon<T>> + Clone,
+ Container: IntoIterator<Item = impl BooleanOps<Scalar = T>> + Clone,
This should work for both Iterator<Item=Polygon> and Iterator<Item=MultiPolygon>, but if you wanted a single iterator, you'd need to do something like Iterator<Item=Box<dyn BooleanOps>>

You can't use impl in trait bounds, so that won't work as written (and it has to be an RTreeObject but that can probably be tagged on.

I think we will actually need separate impls even if I can figure out the above because I'm struggling to imagine how we'd specialise the fold function definition inside unary_union:

it's currently

let fold = |mut accum: MultiPolygon<T>, poly: &Polygon<T>| -> MultiPolygon<T> {
            accum = accum.union(poly);
            accum
        };

but for MultiPolygon it would have to be

let fold = |mut accum: MultiPolygon<T>, mpoly: &MultiPolygon<T>| -> MultiPolygon<T> {
            accum = accum.union(mpoly);
            accum
        };

urschrei · 2024-11-02T00:54:02Z

I wonder if this would work if Polygon and MultiPolygon were enum members and BooleanOps was impl'd on the enum (followed by UnaryUnion).

urschrei · 2024-11-02T16:30:30Z

No point in tying myself in knots playing type Tetris: MultiPolygon doesn't impl RTreeObject anyway (we should probably do that) so it'll have to be a manual op for now (it's not hugely complex: you flat_map the MultiPolygon containers).

I have a more niche use case of limiting the maximum area of unioned polygons, so I'll probably adapt what's here to do that. I don't think baking in support for that level of customization makes sense.

I think this might be tricky to implement either on your side or in geo as the fold operation definitionally has a scalar output: there's no way to skip geometries or end early – those geometries then wouldn't appear in the output at all, which isn't what you want.

urschrei · 2024-11-02T19:35:15Z

For follow-on work, a preliminary parallel version was proposed in georust/rstar#80 (comment):

fn parallel_bottom_up_fold_reduce<T, S, I, F, R>(tree: &RTree<T>, init: I, fold: F, reduce: R) -> S
where
    T: RTreeObject,
    RTreeNode<T>: Send + Sync,
    S: Send,
    I: Fn() -> S + Send + Sync,
    F: Fn(S, &T) -> S + Send + Sync,
    R: Fn(S, S) -> S + Send + Sync,
{
    fn inner<T, S, I, F, R>(parent: &ParentNode<T>, init: &I, fold: &F, reduce: &R) -> S
    where
        T: RTreeObject,
        RTreeNode<T>: Send + Sync,
        S: Send,
        I: Fn() -> S + Send + Sync,
        F: Fn(S, &T) -> S + Send + Sync,
        R: Fn(S, S) -> S + Send + Sync,
    {
        parent
            .children()
            .into_par_iter()
            .fold(init, |accum, child| match child {
                RTreeNode::Leaf(value) => fold(accum, value),
                RTreeNode::Parent(parent) => {
                    let value = inner(parent, init, fold, reduce);

                    reduce(accum, value)
                }
            })
            .reduce(init, reduce)
    }

    inner(tree.root(), &init, &fold, &reduce)
}

geo/src/algorithm/bool_ops/mod.rs

michaelkirk · 2024-11-02T20:41:08Z

type Tetris

The type tetris isn't too bad. But concerningly the test doesn't pass for multipolygon.

I'm not sure if this represents a bug in BooleanOps or in UnaryUnion (or something else): 39a2f62

I'm OK without merging support for MultiPolygons at present, but wanted to highlight this behavior in case it's meaningful to you.

urschrei · 2024-11-02T20:49:46Z

type Tetris

The type tetris isn't too bad. But concerningly the test doesn't pass for multipolygon.

I'm not sure if this represents a bug in BooleanOps or in UnaryUnion (or something else): 39a2f62

I'm OK without merging support for MultiPolygons at present, but wanted to highlight this behavior in case it's meaningful to you.

Oh nice job, I burned a couple of hours on these today. The test result should be the same AFAIK. What's the actual output geometry??

adamreichold · 2024-11-02T21:07:44Z

Some preliminary perf data:

One thing I wonder about the performance is whether the [CachedEnvelope] combinator can help with the R tree construction due the relatively expensive envelope computation?

urschrei · 2024-11-03T00:18:37Z

type Tetris

The type tetris isn't too bad. But concerningly the test doesn't pass for multipolygon.
I'm not sure if this represents a bug in BooleanOps or in UnaryUnion (or something else): 39a2f62
I'm OK without merging support for MultiPolygons at present, but wanted to highlight this behavior in case it's meaningful to you.

Oh nice job, I burned a couple of hours on these today. The test result should be the same AFAIK. What's the actual output geometry??

The correct output should be:

Polygon { exterior:
    LineString([
        Coord { x: 200.38308698623, y: 288.3387931645567 },
        Coord { x: 204.07584924189865, y: 288.2701221168692 },
        Coord { x: 210.99999999939024, y: 291.99999999857494 },
        Coord { x: 209.99999999939024, y: 289.99999999857494 },
        Coord { x: 212.24082540659725, y: 285.47846008694717 },
        Coord { x: 205.630617245422, y: 287.73853614038774 },
        Coord { x: 204.00000000684082, y: 287.0000000060255 },
        Coord { x: 200.38308698623, y: 288.3387931645567 }]), interiors: [] }]

But what we get is:

[Polygon { exterior: 
    LineString([
        Coord { x: 200.38308698977124, y: 288.3387931571061 },
        Coord { x: 204.07584924543988, y: 288.2701221094186 },
        Coord { x: 205.63061724896323, y: 287.73853613293716 },
        Coord { x: 204.00000001038205, y: 287.0000000060255 },
        Coord { x: 200.38308698977124, y: 288.3387931571061 }]), interiors: [] },
Polygon { exterior: LineString([
    Coord { x: 204.0758492379893, y: 288.2701221168692 },
    Coord { x: 210.9999999954809, y: 291.99999999857494 },
    Coord { x: 209.9999999954809, y: 289.99999999857494 },
    Coord { x: 212.24082541013848, y: 285.47846008694717 },
    Coord { x: 212.24082539523732, y: 285.47846009439775 },
    Coord { x: 212.2408254026879, y: 285.47846008694717 },
    Coord { x: 205.63061724896323, y: 287.73853613293716 },
    Coord { x: 205.6306172564138, y: 287.73853614038774 },
    Coord { x: 204.07584926034104, y: 288.2701221094186 },
    Coord { x: 204.07584924543988, y: 288.2701221094186 },
    Coord { x: 204.0758492379893, y: 288.2701221168692 }]), interiors: [] }
];

That's not any kind of correct. Hm.

urschrei · 2024-11-03T01:29:11Z

type Tetris

The type tetris isn't too bad. But concerningly the test doesn't pass for multipolygon.
I'm not sure if this represents a bug in BooleanOps or in UnaryUnion (or something else): 39a2f62
I'm OK without merging support for MultiPolygons at present, but wanted to highlight this behavior in case it's meaningful to you.

Oh nice job, I burned a couple of hours on these today. The test result should be the same AFAIK. What's the actual output geometry??

The correct output should be:

Polygon { exterior:
    LineString([
        Coord { x: 200.38308698623, y: 288.3387931645567 },
        Coord { x: 204.07584924189865, y: 288.2701221168692 },
        Coord { x: 210.99999999939024, y: 291.99999999857494 },
        Coord { x: 209.99999999939024, y: 289.99999999857494 },
        Coord { x: 212.24082540659725, y: 285.47846008694717 },
        Coord { x: 205.630617245422, y: 287.73853614038774 },
        Coord { x: 204.00000000684082, y: 287.0000000060255 },
        Coord { x: 200.38308698623, y: 288.3387931645567 }]), interiors: [] }]

But what we get is:

[Polygon { exterior: 
    LineString([
        Coord { x: 200.38308698977124, y: 288.3387931571061 },
        Coord { x: 204.07584924543988, y: 288.2701221094186 },
        Coord { x: 205.63061724896323, y: 287.73853613293716 },
        Coord { x: 204.00000001038205, y: 287.0000000060255 },
        Coord { x: 200.38308698977124, y: 288.3387931571061 }]), interiors: [] },
Polygon { exterior: LineString([
    Coord { x: 204.0758492379893, y: 288.2701221168692 },
    Coord { x: 210.9999999954809, y: 291.99999999857494 },
    Coord { x: 209.9999999954809, y: 289.99999999857494 },
    Coord { x: 212.24082541013848, y: 285.47846008694717 },
    Coord { x: 212.24082539523732, y: 285.47846009439775 },
    Coord { x: 212.2408254026879, y: 285.47846008694717 },
    Coord { x: 205.63061724896323, y: 287.73853613293716 },
    Coord { x: 205.6306172564138, y: 287.73853614038774 },
    Coord { x: 204.07584926034104, y: 288.2701221094186 },
    Coord { x: 204.07584924543988, y: 288.2701221094186 },
    Coord { x: 204.0758492379893, y: 288.2701221168692 }]), interiors: [] }
];

That's not any kind of correct. Hm.

I think this is a bug in the MultiPolygon union algorithm (not clear if it stems from i_overlay):

Using the test polygons, you should get the same result if you do:

let empty = MultiPolygon::<f64>::new(vec![]);
// multi_poly_1 is [poly1, poly2]
let empty_union = empty.union(&multi_poly_1);

// and
let mp1 = MultiPolygon::new(vec![poly1.clone()]);
let mp2 = MultiPolygon::new(vec![poly2.clone()]);
let mpu = &mp1.union(&mp2);

…but you don't. That's a bug unless I'm missing a detail about how MultiPolygon union is supposed to work.

michaelkirk · 2024-11-03T04:25:05Z

I'm not sure if it's relevant, but it might be interesting that MultiPolygon(vec![poly1, poly2]) is not a valid multipolygon

MULTIPOLYGON(
  ((204 287,206.69670020700084 288.2213844497616,200.38308697914755 288.338793163584,204 287)),
  ((210 290,204.07584923592933 288.2701221108328,212.24082541367974 285.47846008552216,210 290))
)

It'd be nice to be able to handle this case - but it might be a hint towards resolving it.

michaelkirk · 2024-11-03T04:42:15Z

I was able to get the tests passing by making the multipolygon valid: 5900563

Note the usage of assert_relative_eq, some of the output points are slightly shifted (1e-10), but I think that's expected with the grid snapping in i_overlay.

urschrei · 2024-11-03T11:15:34Z

This will pass when #1254 merges.

I was able to get the tests passing by making the multipolygon valid: 5900563

Note the usage of assert_relative_eq, some of the output points are slightly shifted (1e-10), but I think that's expected with the grid snapping in i_overlay.

Urgh I should have checked validity before anything else. Sorry.

…n perf

…e reborrow We now reborrow, but we could also deref: accum = accum.union(&**poly.0);

So we know what the hell we were doing if we ever have to revisit this

michaelkirk reviewed Oct 31, 2024

View reviewed changes

geo/src/algorithm/bool_ops/mod.rs Outdated Show resolved Hide resolved

urschrei changed the title ~~Cascaded union~~ Cascaded / Unary union Oct 31, 2024

urschrei force-pushed the cascaded_union branch from 46d5a17 to 885f428 Compare November 1, 2024 00:14

urschrei force-pushed the cascaded_union branch from 899c782 to 79290f4 Compare November 2, 2024 16:16

urschrei marked this pull request as ready for review November 2, 2024 16:30

urschrei mentioned this pull request Nov 2, 2024

Post-order element traversal georust/rstar#80

Open

adamreichold reviewed Nov 2, 2024

View reviewed changes

geo/src/algorithm/bool_ops/mod.rs Outdated Show resolved Hide resolved

adamreichold reviewed Nov 2, 2024

View reviewed changes

geo/src/algorithm/bool_ops/mod.rs Show resolved Hide resolved

adamreichold reviewed Nov 2, 2024

View reviewed changes

geo/src/algorithm/bool_ops/mod.rs Outdated Show resolved Hide resolved

urschrei mentioned this pull request Nov 3, 2024

Add rstar MultiPolygon compat #1254

Merged

2 tasks

urschrei added 26 commits November 16, 2024 15:19

We're using recursion, probably worth noting

4facd5c

Better top-level doc

7300870

Typo fix

9fa7791

Improve generic names and bounds

c950d8e

Doc improvements

49b0480

Update changelog

ce48cd8

Group constant function parameters to reduce recursion overhead

6b949ae

Switch to CachedEnvelope for tree items to try to improve constructio…

120ab42

…n perf

Implement trait less generally to allows later cloning

4c3e7ac

Switch to boppable collection for unary union impl

34f738c

Fix tests

d35cadd

Typo fix

7465405

Remove need to clone geometries to be unioned before storage in r*-tree

995bd09

Remove deref of immutable expression – the compiler can auto-deref th…

16a7570

…e reborrow We now reborrow, but we could also deref: accum = accum.union(&**poly.0);

Replace manual ObjectRef with rstar ObjectRef

d933095

Implement multithreaded unary union

deb6d5d

Remove redundant Rayon imports

7c0687b

Simplify tree traversal algorithms

0c76949

Ops and driver fns for single-threaded unary union don't need mutability

6388bb4

Document trait bounds in great detail

00106c9

So we know what the hell we were doing if we ever have to revisit this

Multithreading support is now more general than just i_overlay

fd8ae50

Unar union is available for MultiPolygons too

b36ef8c

Clarify what Convert does

dc0ad40

List items don't end with a period

a43a7fb

Note that MultiPolygon can also be stored in an Rtree

31de4dd

Update minimum rstar version to 0.12.2 and update changelogs

83e06bd

urschrei force-pushed the cascaded_union branch from 9f1e643 to 83e06bd Compare November 16, 2024 15:22

urschrei added 2 commits November 17, 2024 15:33

Fix up CHANGES

8eff873

Punctuation

c62ab7a

urschrei mentioned this pull request Nov 18, 2024

union op causes out of bounds error using i_overlay 1.7.4 #1273

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Cascaded / Unary union #1246

Cascaded / Unary union #1246

urschrei commented Oct 31, 2024 •

edited

Loading

lnicola commented Oct 31, 2024 •

edited

Loading

urschrei commented Oct 31, 2024 •

edited

Loading

urschrei commented Nov 1, 2024

dabreegster commented Nov 1, 2024

michaelkirk commented Nov 1, 2024 •

edited

Loading

michaelkirk commented Nov 1, 2024

kylebarron commented Nov 1, 2024

michaelkirk commented Nov 1, 2024 •

edited

Loading

urschrei commented Nov 2, 2024

urschrei commented Nov 2, 2024

urschrei commented Nov 2, 2024

urschrei commented Nov 2, 2024

michaelkirk commented Nov 2, 2024

urschrei commented Nov 2, 2024 •

edited

Loading

adamreichold commented Nov 2, 2024

urschrei commented Nov 3, 2024

urschrei commented Nov 3, 2024

michaelkirk commented Nov 3, 2024 •

edited

Loading

michaelkirk commented Nov 3, 2024 •

edited

Loading

urschrei commented Nov 3, 2024

Cascaded / Unary union #1246

Are you sure you want to change the base?

Cascaded / Unary union #1246

Conversation

urschrei commented Oct 31, 2024 • edited Loading

lnicola commented Oct 31, 2024 • edited Loading

urschrei commented Oct 31, 2024 • edited Loading

urschrei commented Nov 1, 2024

dabreegster commented Nov 1, 2024

michaelkirk commented Nov 1, 2024 • edited Loading

michaelkirk commented Nov 1, 2024

kylebarron commented Nov 1, 2024

michaelkirk commented Nov 1, 2024 • edited Loading

urschrei commented Nov 2, 2024

urschrei commented Nov 2, 2024

urschrei commented Nov 2, 2024

urschrei commented Nov 2, 2024

michaelkirk commented Nov 2, 2024

urschrei commented Nov 2, 2024 • edited Loading

adamreichold commented Nov 2, 2024

urschrei commented Nov 3, 2024

urschrei commented Nov 3, 2024

michaelkirk commented Nov 3, 2024 • edited Loading

michaelkirk commented Nov 3, 2024 • edited Loading

urschrei commented Nov 3, 2024

urschrei commented Oct 31, 2024 •

edited

Loading

lnicola commented Oct 31, 2024 •

edited

Loading

urschrei commented Oct 31, 2024 •

edited

Loading

michaelkirk commented Nov 1, 2024 •

edited

Loading

michaelkirk commented Nov 1, 2024 •

edited

Loading

urschrei commented Nov 2, 2024 •

edited

Loading

michaelkirk commented Nov 3, 2024 •

edited

Loading

michaelkirk commented Nov 3, 2024 •

edited

Loading